library(readxl)
mydata <- read_xlsx("Workload.xlsx")
mydata <- as.data.frame(mydata)
mydata## Training Workload Certification Freq
## 1 1 1 1 37
## 2 1 1 2 15
## 3 1 2 1 52
## 4 1 2 2 18
## 5 2 1 1 14
## 6 2 1 2 41
## 7 2 2 1 21
## 8 2 2 2 46Description of the variables:
# Convert variables to factors with meaningful labels
mydata$Training <- factor(mydata$Training,
levels = c(1, 2),
labels = c("High_TA", "Low_TA"))
mydata$Workload <- factor(mydata$Workload,
levels = c(1, 2),
labels = c("High_W", "Normal_W"))
mydata$Certification <- factor(mydata$Certification,
levels = c(1, 2),
labels = c("Passed", "Failed"))ftable(xtabs(Freq ~ Training + Workload + Certification, data = mydata))## Certification Passed Failed
## Training Workload
## High_TA High_W 37 15
## Normal_W 52 18
## Low_TA High_W 14 41
## Normal_W 21 46model_3way <- glm(Freq ~ Training * Workload * Certification,
family = poisson,
data = mydata)summary(model_3way)##
## Call:
## glm(formula = Freq ~ Training * Workload * Certification, family = poisson,
## data = mydata)
##
## Coefficients:
## Estimate
## (Intercept) 3.61092
## TrainingLow_TA -0.97186
## WorkloadNormal_W 0.34033
## CertificationFailed -0.90287
## TrainingLow_TA:WorkloadNormal_W 0.06514
## TrainingLow_TA:CertificationFailed 1.97738
## WorkloadNormal_W:CertificationFailed -0.15800
## TrainingLow_TA:WorkloadNormal_W:CertificationFailed -0.13239
## Std. Error
## (Intercept) 0.16440
## TrainingLow_TA 0.31378
## WorkloadNormal_W 0.21508
## CertificationFailed 0.30609
## TrainingLow_TA:WorkloadNormal_W 0.40658
## TrainingLow_TA:CertificationFailed 0.43533
## WorkloadNormal_W:CertificationFailed 0.41046
## TrainingLow_TA:WorkloadNormal_W:CertificationFailed 0.57763
## z value Pr(>|z|)
## (Intercept) 21.964 < 2e-16
## TrainingLow_TA -3.097 0.00195
## WorkloadNormal_W 1.582 0.11357
## CertificationFailed -2.950 0.00318
## TrainingLow_TA:WorkloadNormal_W 0.160 0.87271
## TrainingLow_TA:CertificationFailed 4.542 5.57e-06
## WorkloadNormal_W:CertificationFailed -0.385 0.70028
## TrainingLow_TA:WorkloadNormal_W:CertificationFailed -0.229 0.81872
##
## (Intercept) ***
## TrainingLow_TA **
## WorkloadNormal_W
## CertificationFailed **
## TrainingLow_TA:WorkloadNormal_W
## TrainingLow_TA:CertificationFailed ***
## WorkloadNormal_W:CertificationFailed
## TrainingLow_TA:WorkloadNormal_W:CertificationFailed
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 5.4093e+01 on 7 degrees of freedom
## Residual deviance: -7.5495e-15 on 0 degrees of freedom
## AIC: 57.145
##
## Number of Fisher Scoring iterations: 3model_2way <- glm(Freq ~ (Training + Workload + Certification)^2,
family = poisson,
data = mydata)anova(model_3way, model_2way, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: Freq ~ Training * Workload * Certification
## Model 2: Freq ~ (Training + Workload + Certification)^2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 0 0.000000
## 2 1 0.052605 -1 -0.052605 0.8186model_no_TW <- glm(Freq ~ Training + Workload + Certification +
Training:Certification + Workload:Certification,
family = poisson,
data = mydata)
anova(model_2way, model_no_TW, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: Freq ~ (Training + Workload + Certification)^2
## Model 2: Freq ~ Training + Workload + Certification + Training:Certification +
## Workload:Certification
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1 0.052605
## 2 2 0.052607 -1 -1.9451e-06 0.9989model_no_TC <- glm(Freq ~ Training + Workload + Certification +
Training:Workload + Workload:Certification,
family = poisson,
data = mydata)
anova(model_2way, model_no_TC, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: Freq ~ (Training + Workload + Certification)^2
## Model 2: Freq ~ Training + Workload + Certification + Training:Workload +
## Workload:Certification
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1 0.053
## 2 2 49.420 -1 -49.367 2.123e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_no_WC <- glm(Freq ~ Training + Workload + Certification +
Training:Workload + Training:Certification,
family = poisson,
data = mydata)
anova(model_2way, model_no_WC, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: Freq ~ (Training + Workload + Certification)^2
## Model 2: Freq ~ Training + Workload + Certification + Training:Workload +
## Training:Certification
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1 0.05260
## 2 2 0.66261 -1 -0.61 0.4348model_final <- glm(Freq ~ Training + Workload + Certification +
Training:Certification,
family = poisson,
data = mydata)
summary(model_final)##
## Call:
## glm(formula = Freq ~ Training + Workload + Certification + Training:Certification,
## family = poisson, data = mydata)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 3.6643 0.1284 28.541
## TrainingLow_TA -0.9333 0.1995 -4.678
## WorkloadNormal_W 0.2472 0.1290 1.916
## CertificationFailed -0.9921 0.2038 -4.868
## TrainingLow_TA:CertificationFailed 1.9027 0.2857 6.661
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## TrainingLow_TA 2.90e-06 ***
## WorkloadNormal_W 0.0554 .
## CertificationFailed 1.13e-06 ***
## TrainingLow_TA:CertificationFailed 2.73e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 54.09301 on 7 degrees of freedom
## Residual deviance: 0.81243 on 3 degrees of freedom
## AIC: 51.957
##
## Number of Fisher Scoring iterations: 3exp(coef(model_final))## (Intercept) TrainingLow_TA
## 39.0286885 0.3932584
## WorkloadNormal_W CertificationFailed
## 1.2803738 0.3707865
## TrainingLow_TA:CertificationFailed
## 6.7038961tab_TC <- xtabs(Freq ~ Training + Certification, data = mydata)
tab_TC## Certification
## Training Passed Failed
## High_TA 89 33
## Low_TA 35 87OR_TC <- (tab_TC[1,1] * tab_TC[2,2]) / (tab_TC[1,2] * tab_TC[2,1])
OR_TC## [1] 6.703896observed <- mydata$Freq
expected <- fitted(model_final)
# leverage (hat values)
h <- hatvalues(model_final)
# adjusted standardized residuals
residuals_adj <- (observed - expected) / sqrt(expected * (1 - h))
cell_fit <- data.frame(
Training = mydata$Training,
Workload = mydata$Workload,
Certification = mydata$Certification,
Observed = observed,
Expected = round(expected, 3),
Adj_Residual = round(residuals_adj, 3))
cell_fit## Training Workload Certification Observed Expected Adj_Residual
## 1 High_TA High_W Passed 37 39.029 -0.544
## 2 High_TA High_W Failed 15 14.471 0.199
## 3 High_TA Normal_W Passed 52 49.971 0.544
## 4 High_TA Normal_W Failed 18 18.529 -0.199
## 5 Low_TA High_W Passed 14 15.348 -0.496
## 6 Low_TA High_W Failed 41 38.152 0.767
## 7 Low_TA Normal_W Passed 21 19.652 0.496
## 8 Low_TA Normal_W Failed 46 48.848 -0.767LR <- deviance(model_final)
df <- df.residual(model_final)
p_value <- pchisq(LR, df, lower.tail = FALSE)
c(LR = LR, df = df, p_value = p_value)## LR df p_value
## 0.8124270 3.0000000 0.8464924